clear
clc
close all
%% 初始化
n = 1000;%迭代次数
fx = @(x) (x(1)-100)^2+(x(2)-150)^2-10000;%函数
sx2 = @(x1) 300-x1;%边界
syms x1 x2;
x = [x1,x2];
f_sym = fx(x);
dfx_1_1 = diff(diff(f_sym,x1),x1);
dfx_1_2 = diff(diff(f_sym,x1),x2);
dfx_2_2 = diff(diff(f_sym,x2),x2);

dfx = [diff(f_sym,x1);diff(f_sym,x2)];%梯度
H = double([dfx_1_1,dfx_1_2;dfx_1_2,dfx_2_2]);

sample.x = [];
sample.dx = [];
sample.y = [];
x_iteration = repmat(sample,n,1);
%% 绘制计算域
nx = 100;
x_1 = linspace(0,300,nx);
bx2 = sx2(x_1);%边界
x_2 = x_1;
for i = 1:1:nx
    for j = 1:1:nx
        y(i,j) = fx([x_1(i),x_2(j)]);
    end
end
pcolor(x_1,x_2,y);
shading interp;
hold on 
plot(x_1,bx2,'r');

%% 计算
x_iteration(1).x = [150;150];%起始点
x_iteration(1).y = fx(x_iteration(1).x);
x_iteration(1).dx = double(subs(dfx,[x1;x2],x_iteration(1).x));

iteration = [1];
error = norm(x_iteration(1).dx);
Y = x_iteration(1).y;
%figure(2)
%for i = 2:1:n
i = 1;
while (error(end) >= 1e-10) && (i <= n)
    i = i+1;
    % clf
    %梯度下降
    g = x_iteration(i-1).dx;
    alpha = (g'*g)/(g'*H*g);
    x_iteration(i).x = x_iteration(i-1).x - alpha*x_iteration(i-1).dx;

    %边界条件

    x_iteration(i).y = fx(x_iteration(i).x);

    %求梯度
    x_iteration(i).dx = double(subs(dfx,[x1;x2],x_iteration(i).x));

    %绘图
    iteration = [iteration,i];
    error = [error,norm(x_iteration(i).dx)];
    Y = [Y,x_iteration(i).y];
    subplot(2,1,1)
    semilogy(iteration,error);
    title('梯度值');
    subplot(2,1,2)
    semilogy(iteration,Y);
    title('目标函数值');
    % hold on
    % x_all = [x_iteration.x];
    % plot(x_all(1:2:end),x_all(2:2:end),'+','Color','g');
    pause(0.001);
end